CLIA Pancreas_Combined analysis_GL 02082024

library(swimplot) library(coxphf) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(ggplot2) library(scales) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(ComplexHeatmap)

#This analysis is based on CLIA Pancreas Cohort. We will run analysis based on ctDNA and/or CA 19-9 in pre-op, MRD and Surveillance window. MRD window is defined as 2-12 weeks post-surgery and the Surveillance window as 12 weeks post-op if no adjuvant chemotherapy was given and 4 weeks after the end of chemotherapy for patients that received adjuvant chemotherapy.

#ctDNA positivity by stage and window

#Number of Pts at the MRD Window - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.MRD.Window))
circ_datadf <- as.data.frame(circ_data)

total_mrd <- sum(!is.na(circ_data$ctDNA.MRD.Window))
print(total_mrd)

[1] 110

circ_data$ctDNA.MRD.Window <- as.factor(circ_data$ctDNA.MRD.Window)
cont_table_mrd <- table(circ_data$Stage, circ_data$ctDNA.MRD.Window)
print(cont_table_mrd)

     
       0  1
  I   30  8
  II  35  9
  III 14 14

#Number of Pts at the Surveillance Window - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.Surveillance.Window))
circ_datadf <- as.data.frame(circ_data)

total_surv <- sum(!is.na(circ_data$ctDNA.Surveillance.Window))
print(total_surv)

[1] 156

circ_data$ctDNA.Surveillance.Window <- as.factor(circ_data$ctDNA.Surveillance.Window)
cont_table_surv <- table(circ_data$Stage, circ_data$ctDNA.Surveillance.Window)
print(cont_table_surv)

     
       0  1
  I   45 12
  II  42 19
  III 23 15

#Number of Pts with ctDNA anytime post-op - percentage positivity by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.postSx.anytime))
circ_datadf <- as.data.frame(circ_data)

total_anytime <- sum(!is.na(circ_data$ctDNA.postSx.anytime))
print(total_anytime)

[1] 244

circ_data$ctDNA.postSx.anytime <- as.factor(circ_data$ctDNA.postSx.anytime)
cont_table_anytime <- table(circ_data$Stage, circ_data$ctDNA.postSx.anytime)
print(cont_table_anytime)

     
       0  1
  I   49 30
  II  57 45
  III 24 39

#Radiological Recurrence rates by stage and window

#Number of Pts at the MRD Window - percentage of radiological recurrence by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.MRD.Window))
circ_datadf <- as.data.frame(circ_data)

total_mrd <- sum(!is.na(circ_data$Radiological.Recurrence))
print(total_mrd)

[1] 110

circ_data$Radiological.Recurrence <- as.factor(circ_data$Radiological.Recurrence)
cont_table_mrd <- table(circ_data$Stage, circ_data$Radiological.Recurrence)
print(cont_table_mrd)

     
      FALSE TRUE
  I      19   19
  II     28   16
  III    17   11

#Number of Pts at the Surveillance Window - percentage of radiological recurrence by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.Surveillance.Window))
circ_datadf <- as.data.frame(circ_data)

total_surv <- sum(!is.na(circ_data$Radiological.Recurrence))
print(total_surv)

[1] 156

circ_data$Radiological.Recurrence <- as.factor(circ_data$Radiological.Recurrence)
cont_table_surv <- table(circ_data$Stage, circ_data$Radiological.Recurrence)
print(cont_table_surv)

     
      FALSE TRUE
  I      39   18
  II     37   24
  III    21   17

#We run Fisher's exact test for recurrence rates in Stage II vs III in the surveillance window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="I",]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.Surveillance.Window))
circ_datadf <- as.data.frame(circ_data)

circ_data$Stage <- factor(circ_data$Stage, levels = c("II", "III"), labels = c("II", "III"))
circ_data$Radiological.Recurrence <- factor(circ_data$Radiological.Recurrence, levels = c("FALSE", "TRUE"), labels = c("No Recurrence", "Recurrence"))
contingency_table <- table(circ_data$Stage, circ_data$Radiological.Recurrence)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 0.10238, df = 1, p-value = 0.749

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.6764
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.5055338 3.0615479
sample estimates:
odds ratio 
  1.245196 

print(contingency_table)

     
      No Recurrence Recurrence
  II             37         24
  III            21         17

#Number of Pts with ctDNA anytime post-op - percentage of radiological recurrence by stage
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.postSx.anytime))
circ_datadf <- as.data.frame(circ_data)

total_anytime <- sum(!is.na(circ_data$Radiological.Recurrence))
print(total_anytime)

[1] 244

circ_data$Radiological.Recurrence <- as.factor(circ_data$Radiological.Recurrence)
cont_table_anytime <- table(circ_data$Stage, circ_data$Radiological.Recurrence)
print(cont_table_anytime)

     
      FALSE TRUE
  I      43   36
  II     56   46
  III    29   34

#Demographics Table

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data_subset <- circ_data %>%
  select(
    Gender,
    Resectability.Group,
    Surgery,
    NAC,
    Stage,
    Grade,
    MSI,
    KRAS.Detailed,
    BRCA2.Detailed,
    Radiological.Recurrence,
    RecSite,
    DFS_months) %>%
  mutate(
    Gender = factor(Gender, levels = c("FEMALE", "MALE"), labels = c("Female", "Male")),
    Resectability.Group = factor(Resectability.Group, levels = c(1, 2), labels = c("Resectable", "Unresectable")),
    Surgery = factor(Surgery, levels = c(1, 0), labels = c("Yes", "No")),
    NAC = factor(NAC, levels = c("TRUE", "FALSE"), labels = c("Yes", "No")),
    Stage = factor(Stage, levels = c("I", "II", "III", "IV"), labels = c("I", "II", "III", "IV")),
    Grade = factor(Grade, levels = c(1, 2, 3), labels = c("G1", "G2", "G3")),
    MSI = factor(MSI),
    KRAS.Detailed = factor(KRAS.Detailed),
    BRCA2.Detailed = factor(BRCA2.Detailed),
    Radiological.Recurrence = factor(Radiological.Recurrence, levels = c("FALSE", "TRUE"), labels = c("No Recurrence", "Recurrence")),
    RecSite = factor(RecSite),
    DFS_months = as.numeric(DFS_months))
table1 <- circ_data_subset %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)"),
    label = list(
      Gender ~ "Gender",
      Resectability.Group ~ "Resectability",
      Surgery ~ "Surgery Performed",
      NAC ~ "Neoadjuvant Treatment",
      Stage ~ "Overall Stage",
      Grade ~ "Histological Grade",
      MSI ~ "MSI Status",
      KRAS.Detailed ~ "KRAS Status",
      BRCA2.Detailed ~ "BRCA2 Status",
      Radiological.Recurrence ~ "Radiological Recurrence",
      RecSite ~ "Recurrence Site",
      DFS_months ~ "Disease-Free Survival (months)"),) %>%
  bold_labels()
table1

Characteristic	N = 3221
Gender
Female	144 (45%)
Male	178 (55%)
Resectability
Resectable	281 (87%)
Unresectable	41 (13%)
Surgery Performed	271 (84%)
Neoadjuvant Treatment	163 (51%)
Overall Stage
I	87 (27%)
II	119 (37%)
III	73 (23%)
IV	43 (13%)
Histological Grade
G1	25 (12%)
G2	131 (63%)
G3	51 (25%)
Unknown	115
MSI Status
MSI-Hi	1 (0.3%)
MSS	321 (100%)
KRAS Status
	88 (27%)
p.G12D	87 (27%)
p.G12R	46 (14%)
p.G12R, p.A11G	1 (0.3%)
p.G12R, p.G12V	1 (0.3%)
p.G12V	79 (25%)
p.G13D	2 (0.6%)
p.Q61H	13 (4.0%)
p.Q61H, p.Q61L	1 (0.3%)
p.Q61R	3 (0.9%)
p.V14I, p.G12D	1 (0.3%)
BRCA2 Status
	313 (97%)
p.C2212Vfs*17	1 (0.3%)
p.D2983*	1 (0.3%)
p.E1250K	1 (0.3%)
p.E621*	1 (0.3%)
p.G637S, p.V1283Kfs*2	1 (0.3%)
p.K242Q	1 (0.3%)
p.L3172F	1 (0.3%)
p.T2722R	1 (0.3%)
p.T3033Nfs*11	1 (0.3%)
Radiological Recurrence
No Recurrence	148 (46%)
Recurrence	174 (54%)
Recurrence Site
	195 (61%)
Bone	1 (0.3%)
Liver	47 (15%)
Liver, Brain	1 (0.3%)
Liver, Kidney	1 (0.3%)
Liver, Lung	6 (1.9%)
Liver, Peritoneum	1 (0.3%)
Local	31 (9.6%)
Local, Bone	2 (0.6%)
Local, Liver	4 (1.2%)
Local, Lung	6 (1.9%)
Local, Peritoneum	1 (0.3%)
Lung	12 (3.7%)
Lung, Bone	1 (0.3%)
Lymph Node	1 (0.3%)
Mediastinum	1 (0.3%)
Pelvis	1 (0.3%)
Peritoneum	8 (2.5%)
Peritoneum, Bone	1 (0.3%)
Peritoneum, Lung	1 (0.3%)
Disease-Free Survival (months)	14 (0 - 214)
Unknown	51
1 n (%); Median (Range)

fit1 <- as_flex_table(
  table1,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)
fit1

Characteristic	N = 3221
Gender
Female	144 (45%)
Male	178 (55%)
Resectability
Resectable	281 (87%)
Unresectable	41 (13%)
Surgery Performed	271 (84%)
Neoadjuvant Treatment	163 (51%)
Overall Stage
I	87 (27%)
II	119 (37%)
III	73 (23%)
IV	43 (13%)
Histological Grade
G1	25 (12%)
G2	131 (63%)
G3	51 (25%)
Unknown	115
MSI Status
MSI-Hi	1 (0.3%)
MSS	321 (100%)
KRAS Status
	88 (27%)
p.G12D	87 (27%)
p.G12R	46 (14%)
p.G12R, p.A11G	1 (0.3%)
p.G12R, p.G12V	1 (0.3%)
p.G12V	79 (25%)
p.G13D	2 (0.6%)
p.Q61H	13 (4.0%)
p.Q61H, p.Q61L	1 (0.3%)
p.Q61R	3 (0.9%)
p.V14I, p.G12D	1 (0.3%)
BRCA2 Status
	313 (97%)
p.C2212Vfs*17	1 (0.3%)
p.D2983*	1 (0.3%)
p.E1250K	1 (0.3%)
p.E621*	1 (0.3%)
p.G637S, p.V1283Kfs*2	1 (0.3%)
p.K242Q	1 (0.3%)
p.L3172F	1 (0.3%)
p.T2722R	1 (0.3%)
p.T3033Nfs*11	1 (0.3%)
Radiological Recurrence
No Recurrence	148 (46%)
Recurrence	174 (54%)
Recurrence Site
	195 (61%)
Bone	1 (0.3%)
Liver	47 (15%)
Liver, Brain	1 (0.3%)
Liver, Kidney	1 (0.3%)
Liver, Lung	6 (1.9%)
Liver, Peritoneum	1 (0.3%)
Local	31 (9.6%)
Local, Bone	2 (0.6%)
Local, Liver	4 (1.2%)
Local, Lung	6 (1.9%)
Local, Peritoneum	1 (0.3%)
Lung	12 (3.7%)
Lung, Bone	1 (0.3%)
Lymph Node	1 (0.3%)
Mediastinum	1 (0.3%)
Pelvis	1 (0.3%)
Peritoneum	8 (2.5%)
Peritoneum, Bone	1 (0.3%)
Peritoneum, Lung	1 (0.3%)
Disease-Free Survival (months)	14 (0 - 214)
Unknown	51
1n (%); Median (Range)

save_as_docx(fit1, path= "~/Downloads/table1.docx")

#Heatmap with Clinical & Genomics Factors

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data %>% arrange(Stage)
circ_datadf <- as.data.frame(circ_data)

ha <- HeatmapAnnotation(
  Stage = circ_data$Stage,
  Gender = circ_data$Gender,
  Surgery = circ_data$Surgery,
  NAC = circ_data$NAC,
  Grade = circ_data$Grade,
  MSI = circ_data$MSI,
  KRAS = circ_data$KRAS,
  BRCA2 = circ_data$BRCA2,
  ctDNA.MRD.Window = circ_data$ctDNA.MRD.Window,
  ctDNA.Surveillance.Window = circ_data$ctDNA.Surveillance.Window,
  ctDNA.anytime = circ_data$ctDNA.anytime,
  Radiological.Recurrence = circ_data$Radiological.Recurrence,
  
  col = list(Stage = c("I" = "seagreen1", "II" = "orange", "III" = "purple", "IV" = "brown2"),
    Gender = c("FEMALE" = "goldenrod" , "MALE" = "blue4"),
    Surgery = c("0" = "#C1211A", "1" ="#008BCE"),
    NAC = c("FALSE" = "cornflowerblue", "TRUE" ="darkmagenta"),
    Grade = c("0" = "grey", "1" = "coral", "2" ="darkgreen", "3" = "yellow3"),
    MSI = c("MSS" = "grey", "MSI-Hi" ="black"),
    KRAS = c("0" = "grey", "1" ="black"),
    BRCA2 = c("0" = "grey", "1" ="black"),
    ctDNA.MRD.Window = c("1" = "red3", "0" ="blue"),
    ctDNA.Surveillance.Window = c("1" = "red3", "0" ="blue"),
    ctDNA.anytime = c("1" = "red3", "0" ="blue"),
    Radiological.Recurrence = c("TRUE" = "red3", "FALSE" ="blue")
)
)
ht <- Heatmap(matrix(nrow = 0, ncol = length(circ_data$Stage)),show_row_names = FALSE,cluster_rows = F,cluster_columns = FALSE, top_annotation = ha)
pdf("heatmap.pdf",width = 7, height = 7)
draw(ht, annotation_legend_side = "bottom")
dev.off()

null device 
          1 

#Distribution of ctDNA

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!=4,]
circ_data <- circ_data[circ_data$ctDNA.Surg.1<=24,]
circ_datadf <- as.data.frame(circ_data)

hist(circ_data$ctDNA.Surg.1,
     col='gray',
     main='Surgery to first testing ctDNA interval (<6 months)',
     xlab='Weeks from date of surgery to first ctDNA test',
     ylab='ctDNA Samples',
     ylim=c(0,35),
     breaks=seq(0,24,1), xaxp=c(0,24,24))

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 <= 24) %>%
  tally()
print(filtered_data)

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 <= 2) %>%
  tally()
print(filtered_data)

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 >2, ctDNA.Surg.1 <=4) %>%
  tally()
print(filtered_data)

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 >4, ctDNA.Surg.1 <=8) %>%
  tally()
print(filtered_data)

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 >8, ctDNA.Surg.1 <=12) %>%
  tally()
print(filtered_data)

filtered_data <- circ_data %>%
  filter(ctDNA.Surg.1 >12) %>%
  tally()
print(filtered_data)

#DFS by ctDNA at the MRD window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.MRD.Window))
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.Window, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD.Window=0 79     26  33.31   26.15      NA
ctDNA.MRD.Window=1 31     25   5.19    3.38    11.5

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD Window", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     43      19    0.729  0.0535        0.608        0.818
   24     21       3    0.669  0.0594        0.539        0.771

                ctDNA.MRD.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     12.0000       4.0000      22.0000       0.2138       0.0848       0.0779       0.3937 

circ_data$ctDNA.MRD.Window <- factor(circ_data$ctDNA.MRD.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data)

  n= 110, number of events= 51 

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRD.Window1 1.7577    5.7992   0.3021 5.818 5.94e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.Window1     5.799     0.1724     3.208     10.48

Concordance= 0.686  (se = 0.031 )
Likelihood ratio test= 31.88  on 1 df,   p=2e-08
Wald test            = 33.85  on 1 df,   p=6e-09
Score (logrank) test = 42.54  on 1 df,   p=7e-11

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.8 (3.21-10.48); p = 0"

#DFS by ctDNA at the MRD window - Landmark 12 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.MRD.Window))
circ_data$DFS_months=circ_data$DFS_months-3
circ_data <- circ_data[circ_data$DFS_months>=0,]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.Window, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD.Window=0 74     25  30.31   23.15      NA
ctDNA.MRD.Window=1 22     17   4.82    3.37      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD Window", ylab= "Disease Free Survival", xlab="Time from Landmark Timepoint (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     35      20    0.702  0.0568        0.575        0.797
   24     14       3    0.605  0.0723        0.449        0.729

                ctDNA.MRD.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
    12.00000      1.00000     16.00000      0.09626      0.08676      0.00688      0.33271 

circ_data$ctDNA.MRD.Window <- factor(circ_data$ctDNA.MRD.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data)

  n= 96, number of events= 42 

                   coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRD.Window1 1.537     4.652    0.335 4.589 4.46e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.Window1     4.652     0.2149     2.413     8.971

Concordance= 0.639  (se = 0.036 )
Likelihood ratio test= 18.74  on 1 df,   p=1e-05
Wald test            = 21.06  on 1 df,   p=4e-06
Score (logrank) test = 25.26  on 1 df,   p=5e-07

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 4.65 (2.41-8.97); p = 0"

#DFS by ctDNA and NAC combination - four groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60

circ_data$ctDNA.MRD.NAC <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.MRD.NAC = case_when(
    NAC == 1 & ctDNA.MRD.Window == 0 ~ 1,
    NAC == 1 & ctDNA.MRD.Window == 1 ~ 2,
    NAC == 0 & ctDNA.MRD.Window == 0 ~ 3,
    NAC == 0 & ctDNA.MRD.Window == 1 ~ 4
  ))

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.NAC, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.NAC, data = circ_data)

   149 observations deleted due to missingness 
                 n events median 0.95LCL 0.95UCL
ctDNA.MRD.NAC=1 53     18  29.34   26.15      NA
ctDNA.MRD.NAC=2 14     13   5.78    4.63      NA
ctDNA.MRD.NAC=3 26      8     NA   13.14      NA
ctDNA.MRD.NAC=4 17     12   7.92    3.06      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.NAC, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple", "red"), title="DFS - ctDNA MRD Window & NAC", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("NAC & ctDNA(-)", "NAC & ctDNA(+)","Surgery & ctDNA(-)", "Surgery & ctDNA(+)"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.NAC, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

149 observations deleted due to missingness 
                ctDNA.MRD.NAC=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     53       0    1.000  0.0000        1.000        1.000
   24     16      14    0.704  0.0667        0.551        0.813

                ctDNA.MRD.NAC=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           14            0            1            0            1            1 

                ctDNA.MRD.NAC=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     26       0    1.000   0.000        1.000        1.000
   24      5       8    0.542   0.137        0.253        0.762

                ctDNA.MRD.NAC=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           17            0            1            0            1            1 

circ_data$ctDNA.MRD.NAC <- factor(circ_data$ctDNA.MRD.NAC, levels=c("1","2", "3", "4"), labels = c("NAC & ctDNA(-)", "NAC & ctDNA(+)","Surgery & ctDNA(-)", "Surgery & ctDNA(+)"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.NAC, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.NAC, data = circ_data)

  n= 110, number of events= 51 
   (149 observations deleted due to missingness)

                                   coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRD.NACNAC & ctDNA(+)     1.97501   7.20666  0.39548 4.994 5.92e-07 ***
ctDNA.MRD.NACSurgery & ctDNA(-) 0.09668   1.10151  0.42806 0.226    0.821    
ctDNA.MRD.NACSurgery & ctDNA(+) 1.63537   5.13135  0.38792 4.216 2.49e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.NACNAC & ctDNA(+)         7.207     0.1388     3.320    15.645
ctDNA.MRD.NACSurgery & ctDNA(-)     1.102     0.9078     0.476     2.549
ctDNA.MRD.NACSurgery & ctDNA(+)     5.131     0.1949     2.399    10.976

Concordance= 0.692  (se = 0.035 )
Likelihood ratio test= 32.62  on 3 df,   p=4e-07
Wald test            = 34.66  on 3 df,   p=1e-07
Score (logrank) test = 44.13  on 3 df,   p=1e-09

#DFS by ctDNA at the MRD window - Stage analysis

#Analysis at Stage I
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="I",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.Window, data = circ_data)

   44 observations deleted due to missingness 
                    n events median 0.95LCL 0.95UCL
ctDNA.MRD.Window=0 30     12  26.15   13.34      NA
ctDNA.MRD.Window=1  8      7   2.76    1.84      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD Window | Stage I", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(0,24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

44 observations deleted due to missingness 
                ctDNA.MRD.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     30       0    1.000   0.000        1.000        1.000
   24      6      11    0.567   0.107        0.337        0.743

                ctDNA.MRD.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            8            0            1            0            1            1 

circ_data$ctDNA.MRD.Window <- factor(circ_data$ctDNA.MRD.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data)

  n= 38, number of events= 19 
   (44 observations deleted due to missingness)

                     coef exp(coef) se(coef)    z Pr(>|z|)    
ctDNA.MRD.Window1  3.0001   20.0868   0.6608 4.54 5.62e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.Window1     20.09    0.04978     5.501     73.34

Concordance= 0.705  (se = 0.047 )
Likelihood ratio test= 19.95  on 1 df,   p=8e-06
Wald test            = 20.61  on 1 df,   p=6e-06
Score (logrank) test = 35.75  on 1 df,   p=2e-09

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 20.09 (5.5-73.34); p = 0"

#Analysis at Stage II
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="II",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.Window, data = circ_data)

   67 observations deleted due to missingness 
                    n events median 0.95LCL 0.95UCL
ctDNA.MRD.Window=0 35     12  33.31   24.67      NA
ctDNA.MRD.Window=1  9      8   4.63    2.96      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD Window | Stage II", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(0,24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

67 observations deleted due to missingness 
                ctDNA.MRD.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     35       0    1.000   0.000          1.0        1.000
   24     10       9    0.701   0.085          0.5        0.834

                ctDNA.MRD.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            9            0            1            0            1            1 

circ_data$ctDNA.MRD.Window <- factor(circ_data$ctDNA.MRD.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data)

  n= 44, number of events= 20 
   (67 observations deleted due to missingness)

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRD.Window1 2.0295    7.6103   0.5052 4.017 5.89e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.Window1      7.61     0.1314     2.827     20.48

Concordance= 0.686  (se = 0.052 )
Likelihood ratio test= 13.92  on 1 df,   p=2e-04
Wald test            = 16.14  on 1 df,   p=6e-05
Score (logrank) test = 21.69  on 1 df,   p=3e-06

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 7.61 (2.83-20.48); p = 0"

#Analysis at Stage III
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="III",]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.Window, data = circ_data)

   38 observations deleted due to missingness 
                    n events median 0.95LCL 0.95UCL
ctDNA.MRD.Window=0 14      2     NA      NA      NA
ctDNA.MRD.Window=1 14     10   11.4    7.82      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD Window | Stage III", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(0,24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

38 observations deleted due to missingness 
                ctDNA.MRD.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     14       0      1.0   0.000        1.000        1.000
   24      5       2      0.8   0.126        0.409        0.946

                ctDNA.MRD.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           14            0            1            0            1            1 

circ_data$ctDNA.MRD.Window <- factor(circ_data$ctDNA.MRD.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.Window, data = circ_data)

  n= 28, number of events= 12 
   (38 observations deleted due to missingness)

                   coef exp(coef) se(coef)     z Pr(>|z|)   
ctDNA.MRD.Window1 2.153     8.609    0.801 2.688  0.00719 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.Window1     8.609     0.1162     1.791     41.37

Concordance= 0.715  (se = 0.057 )
Likelihood ratio test= 10.08  on 1 df,   p=0.002
Wald test            = 7.22  on 1 df,   p=0.007
Score (logrank) test = 9.73  on 1 df,   p=0.002

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 8.61 (1.79-41.37); p = 0.007"

#DFS by ctDNA at the Surveillance window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Surveillance.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance.Window, data = circ_data)

   103 observations deleted due to missingness 
                              n events median 0.95LCL 0.95UCL
ctDNA.Surveillance.Window=0 110     17     NA      NA      NA
ctDNA.Surveillance.Window=1  46     43   11.4    10.2    14.5

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=12, palette=c("blue","red"), title="DFS - ctDNA Surveillance Window", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

103 observations deleted due to missingness 
                ctDNA.Surveillance.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     87       7    0.934  0.0241        0.866        0.968
   24     45       6    0.853  0.0388        0.757        0.913

                ctDNA.Surveillance.Window=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     21      25    0.457  0.0734       0.3097        0.592
   24      5      16    0.109  0.0459       0.0399        0.217

circ_data$ctDNA.Surveillance.Window <- factor(circ_data$ctDNA.Surveillance.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data)

  n= 156, number of events= 60 
   (103 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Surveillance.Window1  2.386    10.873    0.292 8.172 3.03e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Surveillance.Window1     10.87    0.09197     6.135     19.27

Concordance= 0.763  (se = 0.028 )
Likelihood ratio test= 76.07  on 1 df,   p=<2e-16
Wald test            = 66.78  on 1 df,   p=3e-16
Score (logrank) test = 100.5  on 1 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 10.87 (6.13-19.27); p = 0"

#DFS by ctDNA at the Surveillance window - Stage analysis

#Analysis at Stage I
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="I",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Surveillance.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance.Window, data = circ_data)

   25 observations deleted due to missingness 
                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance.Window=0 45      7     NA      NA      NA
ctDNA.Surveillance.Window=1 12     11   9.21    6.37      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=12, palette=c("blue","red"), title="DFS - ctDNA Surveillance Window | Stage I", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

25 observations deleted due to missingness 
                ctDNA.Surveillance.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     34       4    0.910  0.0431        0.777        0.965
   24     19       1    0.875  0.0538        0.719        0.947

                ctDNA.Surveillance.Window=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      5       7   0.4167  0.1423      0.15247        0.665
   24      1       4   0.0833  0.0798      0.00505        0.311

circ_data$ctDNA.Surveillance.Window <- factor(circ_data$ctDNA.Surveillance.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data)

  n= 57, number of events= 18 
   (25 observations deleted due to missingness)

                              coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Surveillance.Window1  2.4291   11.3489   0.5027 4.832 1.35e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Surveillance.Window1     11.35    0.08811     4.237      30.4

Concordance= 0.751  (se = 0.053 )
Likelihood ratio test= 22.9  on 1 df,   p=2e-06
Wald test            = 23.35  on 1 df,   p=1e-06
Score (logrank) test = 35.49  on 1 df,   p=3e-09

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 11.35 (4.24-30.4); p = 0"

#Analysis at Stage II
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="II",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Surveillance.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance.Window, data = circ_data)

   50 observations deleted due to missingness 
                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance.Window=0 42      7     NA      NA      NA
ctDNA.Surveillance.Window=1 19     17   12.2    10.3    26.7

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=12, palette=c("blue","red"), title="DFS - ctDNA Surveillance Window | Stage II", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

50 observations deleted due to missingness 
                ctDNA.Surveillance.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     32       2    0.949  0.0352        0.811        0.987
   24     14       3    0.846  0.0645        0.665        0.934

                ctDNA.Surveillance.Window=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     10       9    0.526  0.1145       0.2872        0.719
   24      4       6    0.211  0.0935       0.0656        0.410

circ_data$ctDNA.Surveillance.Window <- factor(circ_data$ctDNA.Surveillance.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data)

  n= 61, number of events= 24 
   (50 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Surveillance.Window1 1.9251    6.8557   0.4506 4.272 1.94e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Surveillance.Window1     6.856     0.1459     2.834     16.58

Concordance= 0.732  (se = 0.048 )
Likelihood ratio test= 20.71  on 1 df,   p=5e-06
Wald test            = 18.25  on 1 df,   p=2e-05
Score (logrank) test = 24.47  on 1 df,   p=8e-07

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 6.86 (2.83-16.58); p = 0"

#Analysis at Stage III
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage=="III",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Surveillance.Window, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance.Window, data = circ_data)

   28 observations deleted due to missingness 
                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance.Window=0 23      3     NA      NA      NA
ctDNA.Surveillance.Window=1 15     15   11.4    9.27    15.4

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance.Window, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=12, palette=c("blue","red"), title="DFS - ctDNA Surveillance Window | Stage III", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

28 observations deleted due to missingness 
                ctDNA.Surveillance.Window=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     21       1    0.957  0.0425        0.729        0.994
   24     12       2    0.841  0.0858        0.577        0.947

                ctDNA.Surveillance.Window=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      12.000        6.000        9.000        0.400        0.126        0.165        0.628 

circ_data$ctDNA.Surveillance.Window <- factor(circ_data$ctDNA.Surveillance.Window, levels=c("0","1"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance.Window, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance.Window, data = circ_data)

  n= 38, number of events= 18 
   (28 observations deleted due to missingness)

                              coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Surveillance.Window1  3.0700   21.5415   0.6594 4.656 3.23e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Surveillance.Window1     21.54    0.04642     5.915     78.45

Concordance= 0.809  (se = 0.039 )
Likelihood ratio test= 31.33  on 1 df,   p=2e-08
Wald test            = 21.67  on 1 df,   p=3e-06
Score (logrank) test = 37.41  on 1 df,   p=1e-09

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 21.54 (5.92-78.45); p = 0"

#Multivariate regression for DFS at the Surveillance Window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- subset(circ_data, !is.na(ctDNA.Surveillance.Window))
circ_datadf <- as.data.frame(circ_data)

circ_datadf$Gender <- factor(circ_datadf$Gender, levels = c("FEMALE", "MALE"), labels = c("Female", "Male"))
circ_datadf$AT <- factor(circ_datadf$AT, levels = c("0", "1"), labels = c("No", "Yes"))
circ_datadf$NAC <- factor(circ_datadf$NAC, levels = c("TRUE", "FALSE"), labels = c("Therapy", "Upfront Surgery"))
circ_datadf$KRAS <- factor(circ_datadf$KRAS, levels = c("0", "1"), labels = c("wt", "Mut"))
circ_datadf$Stage <- factor(circ_datadf$Stage, levels = c("I", "II", "III"), labels = c("I", "II", "III"))
circ_datadf$CA.Surveillance.Window <- factor(circ_datadf$CA.Surveillance.Window, levels = c("0", "1"), labels = c("CA19-9 < 37 U/mL", "CA19-9 ≥ 37 U/mL"))
circ_datadf$ctDNA.Surveillance.Window <- factor(circ_datadf$ctDNA.Surveillance.Window, levels=c("0","1"), labels = c("Negative", "Positive"))
surv_object<-Surv(time = circ_datadf$DFS_months, event = circ_datadf$DFS.Event) 
cox_fit <- coxph(surv_object ~ Gender + AT + NAC + KRAS + Stage + CA.Surveillance.Window + ctDNA.Surveillance.Window, data=circ_datadf) 
ggforest(cox_fit, data = circ_datadf, main = "Multivariate Regression Model for DFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#ctDNA Dynamics from MRD window to anytime

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_data <- subset(circ_data, !is.na(ctDNA.Dynamics.anytime))
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Dynamics.anytime, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics.anytime, data = circ_data)

                          n events median 0.95LCL 0.95UCL
ctDNA.Dynamics.anytime=1 48      5     NA      NA      NA
ctDNA.Dynamics.anytime=2  3      2  14.72    1.15      NA
ctDNA.Dynamics.anytime=3 19     18   8.02    6.21    24.7
ctDNA.Dynamics.anytime=4 15     14   6.67    3.84      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics.anytime, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple", "red"), title="DFS - ctDNA Dynamics MRD to anytime", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("Persistently negative", "Converted negative","Converted positive", "Persistently positive"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics.anytime, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics.anytime=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     48       0    1.000  0.0000        1.000        1.000
   24     17       4    0.906  0.0451        0.767        0.964

                ctDNA.Dynamics.anytime=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            3            0            1            0            1            1 

                ctDNA.Dynamics.anytime=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     19       0    1.000  0.0000       1.0000         1.00
   24      4      15    0.211  0.0935       0.0656         0.41

                ctDNA.Dynamics.anytime=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           15            0            1            0            1            1 

circ_data$ctDNA.Dynamics.anytime <- factor(circ_data$ctDNA.Dynamics.anytime, levels=c("1", "2", "3", "4"), labels = c("Persistently negative", "Converted negative", "Converted positive", "Persistently positive"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics.anytime, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics.anytime, data = circ_data)

  n= 85, number of events= 39 

                                               coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Dynamics.anytimeConverted negative     2.3940   10.9569   0.8461 2.829  0.00466 ** 
ctDNA.Dynamics.anytimeConverted positive     2.6610   14.3099   0.5076 5.242 1.58e-07 ***
ctDNA.Dynamics.anytimePersistently positive  3.3164   27.5597   0.5471 6.062 1.35e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                            exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Dynamics.anytimeConverted negative        10.96    0.09127     2.087     57.53
ctDNA.Dynamics.anytimeConverted positive        14.31    0.06988     5.292     38.70
ctDNA.Dynamics.anytimePersistently positive     27.56    0.03628     9.432     80.53

Concordance= 0.802  (se = 0.034 )
Likelihood ratio test= 59.42  on 3 df,   p=8e-13
Wald test            = 38.02  on 3 df,   p=3e-08
Score (logrank) test = 67.32  on 3 df,   p=2e-14

#ctDNA Dynamics from MRD window to anytime - Landmark 12 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data$DFS_months=circ_data$DFS_months-3
circ_data <- circ_data[circ_data$DFS_months>=0,]
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
circ_data <- subset(circ_data, !is.na(ctDNA.Dynamics.anytime))
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.Dynamics.anytime, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics.anytime, data = circ_data)

                          n events median 0.95LCL 0.95UCL
ctDNA.Dynamics.anytime=1 48      5     NA      NA      NA
ctDNA.Dynamics.anytime=2  2      1  11.72      NA      NA
ctDNA.Dynamics.anytime=3 19     18   5.02    3.21    21.7
ctDNA.Dynamics.anytime=4 13     12   4.69    2.19      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics.anytime, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple", "red"), title="DFS - ctDNA Dynamics MRD to anytime", ylab= "Disease Free Survival", xlab="Time From Landmark analysis (Months)", legend.labs=c("Persistently negative", "Converted negative","Converted positive", "Persistently positive"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics.anytime, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics.anytime=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     48       0    1.000  0.0000        1.000        1.000
   24     12       4    0.906  0.0451        0.767        0.964

                ctDNA.Dynamics.anytime=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            2            0            1            0            1            1 

                ctDNA.Dynamics.anytime=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     19       0    1.000  0.0000       1.0000        1.000
   24      2      17    0.105  0.0704       0.0178        0.284

                ctDNA.Dynamics.anytime=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           13            0            1            0            1            1 

circ_data$ctDNA.Dynamics.anytime <- factor(circ_data$ctDNA.Dynamics.anytime, levels=c("1", "2", "3", "4"), labels = c("Persistently negative", "Converted negative", "Converted positive", "Persistently positive"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics.anytime, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics.anytime, data = circ_data)

  n= 82, number of events= 36 

                                               coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Dynamics.anytimeConverted negative     1.6947    5.4452   1.1025 1.537    0.124    
ctDNA.Dynamics.anytimeConverted positive     2.6836   14.6381   0.5075 5.288 1.24e-07 ***
ctDNA.Dynamics.anytimePersistently positive  3.2291   25.2581   0.5580 5.787 7.15e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                            exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Dynamics.anytimeConverted negative        5.445    0.18365    0.6274     47.26
ctDNA.Dynamics.anytimeConverted positive       14.638    0.06831    5.4138     39.58
ctDNA.Dynamics.anytimePersistently positive    25.258    0.03959    8.4618     75.39

Concordance= 0.802  (se = 0.035 )
Likelihood ratio test= 55.12  on 3 df,   p=6e-12
Wald test            = 36.56  on 3 df,   p=6e-08
Score (logrank) test = 64.09  on 3 df,   p=8e-14

#DFS by ctDNA and ACT combination in NAC-treated pts - four groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- circ_data[circ_data$NAC=="TRUE",]

circ_data$ctDNA.MRD.ACT <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.MRD.ACT = case_when(
    AT == 1 & ctDNA.MRD.Window == 0 ~ 1,
    AT == 1 & ctDNA.MRD.Window == 1 ~ 3,
    AT == 0 & ctDNA.MRD.Window == 0 ~ 2,
    AT == 0 & ctDNA.MRD.Window == 1 ~ 4
  ))

circ_data <- subset(circ_data, !is.na(ctDNA.MRD.ACT))
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.ACT, data = circ_data)

                 n events median 0.95LCL 0.95UCL
ctDNA.MRD.ACT=1 36     17  29.34   10.91      NA
ctDNA.MRD.ACT=2 17      1     NA      NA      NA
ctDNA.MRD.ACT=3 12     12   5.78    4.63      NA
ctDNA.MRD.ACT=4  2      1   2.96    2.96      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=12, palette=c("blue","green","orange", "red"), title="DFS - Neoadjuvant Setting w/o Adjuvant Treatment", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("NAC/ctDNA(-)/ACT", "NAC/ctDNA(-)/No ACT", "NAC/ctDNA(+)/ACT", "NAC/ctDNA(+)/No ACT"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.ACT, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD.ACT=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     36       0    1.000  0.0000        1.000        1.000
   24     11      13    0.618  0.0833        0.434        0.757

                ctDNA.MRD.ACT=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     17       0    1.000  0.0000        1.000        1.000
   24      5       1    0.917  0.0798        0.539        0.988

                ctDNA.MRD.ACT=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           12            0            1            0            1            1 

                ctDNA.MRD.ACT=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            2            0            1            0            1            1 

circ_data$ctDNA.MRD.ACT <- factor(circ_data$ctDNA.MRD.ACT, levels=c("2","1","3","4"), labels = c("NAC/ctDNA(-)/No ACT", "NAC/ctDNA(-)/ACT", "NAC/ctDNA(+)/ACT", "NAC/ctDNA(+)/No ACT"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD.ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD.ACT, data = circ_data)

  n= 67, number of events= 31 

                                   coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRD.ACTNAC/ctDNA(-)/ACT     1.989     7.305    1.034 1.923  0.05450 .  
ctDNA.MRD.ACTNAC/ctDNA(+)/ACT     3.680    39.658    1.052 3.497  0.00047 ***
ctDNA.MRD.ACTNAC/ctDNA(+)/No ACT  3.164    23.672    1.428 2.216  0.02667 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                 exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRD.ACTNAC/ctDNA(-)/ACT        7.305    0.13688    0.9623     55.46
ctDNA.MRD.ACTNAC/ctDNA(+)/ACT       39.658    0.02522    5.0422    311.92
ctDNA.MRD.ACTNAC/ctDNA(+)/No ACT    23.672    0.04224    1.4421    388.58

Concordance= 0.744  (se = 0.033 )
Likelihood ratio test= 28.04  on 3 df,   p=4e-06
Wald test            = 23.04  on 3 df,   p=4e-05
Score (logrank) test = 35.54  on 3 df,   p=9e-08

#DFS by ctDNA and ACT combination in upfront surgery pts - four groups

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]
circ_data <- circ_data[circ_data$NAC=="FALSE",]

circ_data$ctDNA.MRD.ACT <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.MRD.ACT = case_when(
    AT == 1 & ctDNA.MRD.Window == 0 ~ 1,
    AT == 1 & ctDNA.MRD.Window == 1 ~ 3,
    AT == 0 & ctDNA.MRD.Window == 0 ~ 2,
    AT == 0 & ctDNA.MRD.Window == 1 ~ 4
  ))

circ_data <- subset(circ_data, !is.na(ctDNA.MRD.ACT))
circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~ctDNA.MRD.ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD.ACT, data = circ_data)

                 n events median 0.95LCL 0.95UCL
ctDNA.MRD.ACT=1 22      8     NA   13.14      NA
ctDNA.MRD.ACT=2  4      0     NA      NA      NA
ctDNA.MRD.ACT=3 13      9   7.92    3.06      NA
ctDNA.MRD.ACT=4  4      3   8.51    1.84      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD.ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","orange", "red"), title="DFS - Neoadjuvant Setting w/o Adjuvant Treatment", ylab= "Disease Free Survival", xlab="Months from surgery", legend.labs=c("Surgery/ctDNA(-)/ACT", "Surgery/ctDNA(-)/No ACT", "Surgery/ctDNA(+)/ACT", "Surgery/ctDNA(+)/No ACT"), legend.title="")

summary(KM_curve, times= c(0, 24))

Call: survfit(formula = surv_object ~ ctDNA.MRD.ACT, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD.ACT=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     22       0    1.000   0.000         1.00        1.000
   24      5       8    0.509   0.136         0.23        0.733

                ctDNA.MRD.ACT=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            4            0            1            0            1            1 

                ctDNA.MRD.ACT=3 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0           13            0            1            0            1            1 

                ctDNA.MRD.ACT=4 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
           0            4            0            1            0            1            1 

circ_data$ctDNA.MRD.ACT <- factor(circ_data$ctDNA.MRD.ACT, levels=c("2","1","3","4"), labels = c("Surgery/ctDNA(-)/No ACT", "Surgery/ctDNA(-)/ACT", "Surgery/ctDNA(+)/ACT", "Surgery/ctDNA(+)/No ACT"))
cox_fit <- coxphf(surv_object ~ ctDNA.MRD.ACT, data=circ_data) 

Warning: Convergence in estimating profile likelihood CI or p-values not attained for all variables.
Consider re-run with smaller maxstep and larger maxit.

summary(cox_fit)

coxphf(formula = surv_object ~ ctDNA.MRD.ACT, data = circ_data)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                                          coef se(coef) exp(coef) lower 0.95 upper 0.95     Chisq          p
ctDNA.MRD.ACTSurgery/ctDNA(-)/ACT    0.7995539 1.548934  2.224548  0.2672538   289.6627 0.3766586 0.53939693
ctDNA.MRD.ACTSurgery/ctDNA(+)/ACT    1.9954580 1.539318  7.355571  0.9185990   951.4111 3.4497120 0.06326279
ctDNA.MRD.ACTSurgery/ctDNA(+)/No ACT 2.0956260 1.604608  8.130529  0.7739581  1098.1239        NA         NA

Likelihood ratio test=9.129411 on 3 df, p=0.02761887, n=43
Wald test = 7.832561 on 3 df, p = 0.04960194

Covariance-Matrix:
                                     ctDNA.MRD.ACTSurgery/ctDNA(-)/ACT ctDNA.MRD.ACTSurgery/ctDNA(+)/ACT ctDNA.MRD.ACTSurgery/ctDNA(+)/No ACT
ctDNA.MRD.ACTSurgery/ctDNA(-)/ACT                             2.399197                          2.258354                             2.255400
ctDNA.MRD.ACTSurgery/ctDNA(+)/ACT                             2.258354                          2.369501                             2.251362
ctDNA.MRD.ACTSurgery/ctDNA(+)/No ACT                          2.255400                          2.251362                             2.574766

#DFS by KRAS mutations

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("CLIA Pancreas_GL.csv")
circ_data <- circ_data[circ_data$Include.analysis==TRUE,]
circ_data <- circ_data[circ_data$Surgery==1,]
circ_data <- circ_data[circ_data$Stage!="IV",]

circ_data$KRAS.status <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(KRAS.status = case_when(
    KRAS.Detailed == "" ~ "1",
    KRAS.Detailed == "p.G12D" ~ "2",
    KRAS.Detailed == "p.G12V" ~ "3",
    KRAS.Detailed == "p.G12R" ~ "4"
  ))

circ_data <- subset(circ_data, !is.na(KRAS.status))
tally <- table(circ_data$KRAS.status)
print(tally)

 1  2  3  4 
73 71 64 36 

circ_data$DFS_months[circ_data$DFS_months > 60] <- 60
survfit(Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)~KRAS.status, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event) ~ 
    KRAS.status, data = circ_data)

               n events median 0.95LCL 0.95UCL
KRAS.status=1 73     28   33.8    24.7      NA
KRAS.status=2 71     35   21.2    15.6      NA
KRAS.status=3 64     43   15.2    11.5    22.2
KRAS.status=4 36     16   31.3    17.5      NA

surv_object <-Surv(time = circ_data$DFS_months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ KRAS.status, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue", "green", "purple", "red"), title="DFS - KRAS mutations", ylab= "Disease Free Survival", xlab="Months from Surgery", legend.labs=c("wt","G12D", "G12V", "G12R"), legend.title="")

summary(KM_curve, times= c(0, 12, 24))

Call: survfit(formula = surv_object ~ KRAS.status, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                KRAS.status=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     73       0    1.000  0.0000        1.000        1.000
   12     47      17    0.749  0.0529        0.627        0.836
   24     28       6    0.635  0.0624        0.500        0.743

                KRAS.status=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     71       0    1.000  0.0000        1.000        1.000
   12     36      24    0.645  0.0587        0.518        0.747
   24     16       8    0.468  0.0692        0.329        0.596

                KRAS.status=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     64       0    1.000  0.0000        1.000        1.000
   12     31      24    0.595  0.0645        0.457        0.708
   24     13      13    0.323  0.0666        0.199        0.454

                KRAS.status=4 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0     36       0    1.000  0.0000        1.000        1.000
   12     24       7    0.793  0.0700        0.613        0.896
   24     14       5    0.588  0.0947        0.382        0.746

circ_data$KRAS.status <- factor(circ_data$KRAS.status, levels =c("1", "2", "3", "4"), labels =c("wt","G12D", "G12V", "G12R"))
cox_fit <- coxph(surv_object ~ KRAS.status, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ KRAS.status, data = circ_data)

  n= 244, number of events= 122 

                   coef exp(coef) se(coef)     z Pr(>|z|)   
KRAS.statusG12D 0.52284   1.68682  0.25487 2.051  0.04023 * 
KRAS.statusG12V 0.76551   2.15009  0.24358 3.143  0.00167 **
KRAS.statusG12R 0.09868   1.10372  0.31369 0.315  0.75308   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
KRAS.statusG12D     1.687     0.5928    1.0236     2.780
KRAS.statusG12V     2.150     0.4651    1.3339     3.466
KRAS.statusG12R     1.104     0.9060    0.5968     2.041

Concordance= 0.58  (se = 0.026 )
Likelihood ratio test= 12.22  on 3 df,   p=0.007
Wald test            = 11.98  on 3 df,   p=0.007
Score (logrank) test = 12.42  on 3 df,   p=0.006